MICROCOPY  RLSOIUIION  TLM  CHARI 

NAMONAi  HliRfAU  l»«  • ' ANI  ■•>•'! 


February  1977 
J77-75008-TR-1 


A BALANCED  EXPANSION  TECHNIQUE 


FOR  CONSTRUCTING  ACCURATE 


FINITE  DIFFERENCE  ADVECTION  SCHEMES* 


Robert  K.-C.  Chan 


Approved  for  Public  Release 
Distribution  Unlimited 


This  research  was  sponsored  under 
Contract  N00014-76-C-0455 
Office  of  Naval  Research 


*Submitted  to  The  Journal  of  Computational  Physics  for  publication 


ABSTRACT 


A new  procedure,  called  the  Balanced  Expansion  Technique 
(BET),  is  employed  to  construct  accurate  finite  difference 
advection  schemes  that,  for  the  model  equation  considered,  are 
neutrally  stable.  By  applying  BET  systematically,  the  phase 
error  can  be  made  as  small  as  one  wishes.  Test  calculations 
with  one-dimensional  problems  have  confirmed  the  expected 
accuracy  of  these  methods. 
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I . INTRODUCTION 

Most  finite  difference  methods  for  solving  the  advection 
equation  (Eq . (1))  suffer  from  some  or  all  of  the  following 
difficulties:  (1)  amplitude  error,  (2)  phase  error,  and 

(3)  Courant  condition.  Methods  for  simulating  the  advection- 
dif fusion  process  (see  Eq.  (26))  may,  in  addition,  be  plagued 
by  the  cell  Reynolds  number  limit  {1}. 

For  pure  convection  at  a constant  velocity,  it  is  well- 
known  that  all  Fourier  component  solutions  travel  at  the  same 
speed  (i.e.,  nondispersive)  and  their  amplitudes  remain 
unchanged  at  all  times.  Many  finite  difference  schemes,  how- 
ever, produce  significant  amount  of  errors  in  both  the  amplitude 
and  phase  of  component  solutions.  Some  methods  are  free  from 
amplitude  errors,  but  do  not  preserve  the  phase  well.  While 
unstable  schemes  that  systematically  magnify  the  amplitudes 
are  obviously  useless,  those  deriving  their  numerical  stability 
from  nonphysical  damping,  which  overwhelms  other  important 
terms  in  the  equation,  are  also  rather  limited  in  their  range 
of  applications.  Accuracy  in  phase  preservation  for  those 
components  in  the  physically  important  range  of  wavelengths  is 
equally  important.  Numerical  dispersion  of  such  components 
can  seriously  degrade  the  quality  of  the  solution,  and,  like 
the  amplitude  error,  though  less  obviously,  is  a major  source  of 
inaccuracy. 


Besides  amplitude  and  phase  distortions,  explicit 
methods  are  usually  limited  by  the  Courant  condition  which 
restricts  the  maximum  allowable  time  step  such  that  no  fluid 
particle  travels  a distance  greater  than  one  mesh  spacing 
within  a time  increment.  Some  implicit  schemes  are  stable  for 
arbitrarily  large  time  steps,  but  at  the  expense  of  accuracy. 
Thus,  with  these  difficulties,  accurate  computation  of  flows 
at  high  Reynolds  number  appears  to  be  a hopeless  task  in 
computational  fluid  dynamics.  The  reader  is  referred  to 
Roache  {1}  for  more  detailed  discussion  on  the  various  existing 
finite  difference  schemes. 

This  paper  presents  a fundamentally  different  procedure 
for  developing  a set  of  relatively  simple  finite  difference 
advection  schemes  that  lead  to  zero  amplitude  error,  while  the 
phase  distortion  can  be  made  as  small  as  one  wishes.  Using 
this  procedure,  called  the  Balanced  Expansion  Technique  (BET), 
four  different  advection  schemes  have  been  constructed  as 
examples.  These  schemes  can  be  generalized  to  allow  the  use 
of  arbitrarily  large  time  steps,  without  sacrificing  accuracy, 
and  to  include  diffusion. 

In  the  sections  that  follow,  the  derivation  of  each 
scheme  is  given . To  demonstrate  the  usefulness  of  the  new 
schemes,  they  are  tested  with  transient  and  steady-state 
problems  with  known  solutions. 
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II.  THE  BALANCED  EXPANSION  TECHNIQUE 


\ 

\ 

The  basic  BET  procedure  is  best  illustrated  with  a few 
examples.  The  model  equation  that  concerns  us  is 

\ 

where  $ is  a passive  scalar  quantity,  x is  the  one-dimensional 
space  coordinate,  and  t is  the  time.  The  constant  advection 
velocity  u is  taken  to  be  positive  without  loss  of  generality. 

x V 

Throughout  this  paper*,  Eq.  (1)  will  be  interpreted  in  terms 
of  fluid  flows  although  it  plays  an  equally  important  role  in 
many  other  branches  of  physical  sciences. 

To  discretize  Eq.  (1),  a one-dimensional  finite  differ- 
ence mesh  with  constant  spacing  fix  will  be  used.  The  subscript 
j and  superscript  n are  so  defined  that  represents  the 
value  of  ♦ at  x ■ Jfix  and  t » nfit,  where  fit  is  the  constant 
time  increment.  In  what  follows,  several  finite  difference 
advection  schemes  will  be  derived,  using  the  general  BET 
procedure . 

2.1  The  A1  Scheme 

The  objective  of  any  discrete  representation  of  Eq.  (1) 
may  be  stated  as:  Given  the  values  of  (for  all  j),  at  the 

time  level  n and  perhaps  several  previous  time  levels,  find 
an  approximation  for  the  new  values  (for  all  J).  This 
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objective  is  usually  accomplished  by  writing  finite  difference 
analogs  for  the  partial  derivatives,  i.e.,  3 4>/ 3 1 and  3(J>/3x, 
in  Eq.  (1).  Recently,  some  very  high  order  schemes  have  been 
developed  for  evaluating  these  derivatives  {1} . In  BET,  a 
drastically  different  approach  is  taken.  No  attempt  is  made 
to  approximate  the  partial  derivatives.  Rather,  one  just 
uses  the  physical  meaning  behind  Eq.  (1)  and  then  proceeds 
to  obtain  an  algebraic  expression  for  through  careful 

mathematical  manipulations. 

As  is  well-known,  Eq.  (1)  implies  the  preservation  of 
the  value  of  4 along  the  characteristic  line  dx/dt  = u in 
the  (x,t)  space.  The  general  solution  is  $(x,t)  = F(x  - ut), 
which  leads  to  the  following  relations 

♦(x,t  + 6t)  * t(x  - u6t,t)  (2) 

♦ (x  + u6t,t)  - $(x,t  - 6t)  . (3) 

For  convenience,  we  introduce  the  Courant  number  o as 


„ - 

“ * "sr 


(4) 


Referring  to  Fig.  1,  Eq.  (2)  states  that  equals  the 

value  of  4>n  located  at  a distance  u6t  (or  o6x)  to  the  left 
of  $ ® . Similarly,  by  Eq.  (3),  is  numerically  the  same 

as  $n  situated  at  a distance  afix  to  the  right  of  $ ® . 

As  shown  in  Fig.  1,  the  positions  of  the  six  quantities 
^j-1’  *j-l’  ♦jli*  j . and  <t> are  symmetric  about  that 

of  $0,  which  is  at  the  midpoint  between  j and  j-1.  Making 
Taylor  series  expansions  about  we  have 
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,n-l 


♦o  + 


’l 

2 + aj  fix  ♦x  + 2 , 


fix2  <t> 


MW 

JT  (5  + “f  8x3  *XXX  * TT  (5  + «f 

A (i  ♦ «f 8x5  *<5) 


XX 


fix4  <J> 


xxxx 


(5) 


,n+l  . 

Vi " ♦© 


8x  *x  + TT  I * “f  81,2  ♦; 


■[**“) 

- TT  (3  + af  8x3  *xxx  + A (5  + °f  8x4  ♦ 

- rr  (I + “f 8x5  ♦“> 


xxxx 


(6) 


In  the  equations  above,  the  subscript  x means  partial  deriva- 
tives of  $,  evaluated  at  the  same  point  as  <f>Q,  and  is 

a shorthand  notation  for  the  fifth  derivative  4> 


XXXXX • 


Subtracting  Eq.  (6)  from  Eq.  (5),  we  obtain' 

♦j’1  • *”-i  ■ 2(s  + °)  8x  *x  * Ti  (5  + °f  8x3  ♦ 

♦AfW*5  *<5)  + - • • 


XXX 


(7) 


Similarly,  by  containing  expansions  of  and  about  ♦Q, 

we  have 

♦3  ~ ♦S-i  * 2(l)6x  *x  + rr  (if 6x3  ♦xxx  + rr  (if 6x5  ♦<5) 


♦ . . 


(8) 
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The  A1  scheme  involves  only  two  space  points  (j  and 
j-1),  but  three  time  levels.  Thus  it  can  not  be  used  at  the 
first  step  of  time  integration.  The  two-level  A3  scheme,  to 
be  described  shortly,  may  be  used  as  a starter  for  this  scheme. 
Although  appears  on  the  right  side  of  Eq . (10),  the  A1 

scheme  is  a fully  explicit  method  if  is  computed  sequen- 

tially in  the  downstream  direction.  The  radius  of  convergence 
for  the  infinite  series  on  the  right  side  of  Eq.  (10)  is 

0 < a < 1 . (11) 

Thus,  for  Eq.  ^10)  to  be  valid,  we  are  limited  by  the  usual 
Courant  condition.  This  restriction,  however,  can  be  removed 
in  a generalized  A1  scheme  as  will  be  presented  later. 


2.2 


The  A2  Scheme 


The  six  quantities  (|>”_2,  <fr“_lf  ^+1 . ♦*},  and  <|>“+1 


are  used  to  construct  the  A2  scheme  (Fig.  2).  Expanding  these 


variables  in  Taylor  series  about  we  can  form  the  three 


symmetric  (or  balanced)  pairs 
tn 


.n 


♦3+1  - *U  - 2 


i ax  ♦x  + 


nr 


6x2  4>  + Jr 

yxxx  5T 


(if  ♦<•> 


(12) 


♦“ " ♦”-! " 2(i)  6x  *x  + ir  (I)  6x3  ♦xxx + rr  (if 6x5  ♦(5) 


(13) 


*J+1  ■ ♦j-1  " 2(i  “ a)6x  *x  + tt  (i  " a)  6x3  * 


xxx 


Again,  by  eliminating  <J>  and  6 among  these  three  equations, 

A XX  X 

we  have 


♦S+1  ■ ♦ b4^-2  - *U  ♦ ■>,[♦;  - *;.i) 

+ a(i  - a2) ( 1 - 2a) ( 2 - a)6x5  <fr(5) 


where 

bj  = | o(l  - a)(  1 - 2a) 

b2  = | d - 2a)(l  + a) ( 2 - a)  . 


(15) 


If  the  truncation  errors  are  dropped,  Eq . (15)  is  called  the 
A2  scheme . 

Similar  to  Al,  this  formula  involves  three  time  levels 
and  is  completely  explicit.  Again,  for  Eq . (15)  to  be  valid, 
the  Courant  condition  Eq . (11)  must  be  satisfied. 


2 . 3 The  A3  Scheme 

This  scheme  is  constructed  by  expanding  a different 
set  of  six  quantities,  i.e.,  ♦j-i*  4>j+1*  ♦”»  and 

about  (Fig.  3)  and  forming  the  balanced  pairs 
- 4>j _2 * " ^j-i»  and  4>j+1  “ $j_i*  Upon  elimination  of 

$x  and  $xxx  among  the  three  resulting  equations,  we  have 

.n+1  « An  j.  „ Ln  a“4’1)  a „ (a"  An+l\ 

♦j  *J-1  + Cl[*3  ~ ♦j-lJ  + C2  (^j-2  " *j+lj 

+ 3I5  a(!  - a)(2  - a)«x5  ♦(5)  + . . . (16) 
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tm&grur.r 


where 


c = 2 - « 

C1  2(1  + ct) 


a 

!(3  - a] 


By  dropping  the  truncation  errors,  Eq . (16)  is  the  A3  scheme. 
This  method  is  implicit  but  involves  only  two  time  levels, 
which  is  very  useful  for  starting  other  multi-level  schemes 
such  as  A1  and  A2 . For  the  A3  scheme  we  also  require 
0 * a < 1. 


2 .4  The  Higher  Order  A4  Scheme 

From  the  principle  illustrated  so  far,  it  is  apparent 
that  higher  order  schemes  can  be  constructed  by  involving 
more  balanced  pairs  of  discretized  <J>  values,  with  the  center 
of  each  pair  located  at  the  same  point  <PQ . One  possibility 
is  to  combine  A1  and  A2 . The  two  pairs  (f1?  - <$>”_,  and 
4>j  - 4>j  are  shared  by  both  schemes.  In  fact  Eqs . (8) 

and  (9)  are  identical  with  Eqs.  (13)  and  (14),  respectively. 
Thus,  by  combining  Eqs.  (7)-(9)  and  Eqs.  (12)-(14),  we  have 
four  linearly  independent  equations  from  which  we  obtain 

♦r  - + - ♦;-*) 


to  - *;-i] 


+ 0(6x') 


(17) 


where 


i 


rt  = a2(l  - 2a) 

2 3(2  + a) 

d3  - (2  - a) ( 1 - 2a) 

Dropping  the  truncation  errors  terms,  Eq . (17)  is  the  very 

accurate  A4  scheme  for  calculating  . Like  the  A1  scheme, 

this  method  involves  three  time  levels,  but  is  explicit  if 

is  computed  in  the  direction  of  increasing  value  of  j. 

From  the  development  above,  one  could  use  one  more 

balanced  pair  of  <J>  values  to  obtain  an  expression  for  4>j+1 

9 

which  has  truncation  errors  of  order  6x  , the  new  pair  being 
$"“2  and  and  the  process  can  be  systematically  extended 

to  any  order  we  wish.  The  more  lower-order  odd  derivatives 
of  <J>  that  are  eliminated  in  the  BET  procedure,  the  less  phase 
error  will  be  introduced  by  the  scheme.  As  will  become  evident 
in  the  following  section,  the  A4  scheme  is  adequate  for  most 
applications . 


13 


III.  STABILITY  AND  PHASE  ERRORS 


By  von  Neumann  stability  analysis  {1},  all  the  schemes 
presented  in  the  preceding  section  are  found  to  be  neutrally 
stable  for  Fourier  component  solutions  of  all  wavelengths. 

That  is,  the  amplitudes  of  the  component  solutions  are  exactly 
preserved  by  these  schemes.  This  desirable  property  is  a 
direct  consequence  of  employing  the  balanced  expansion  pro- 
cedure in  which  the  damping  terms  (i.e.,  all  the  even  deriva- 
tives of  <f>  in  the  Taylor  series)  are  cancelled  identically. 
Thus,  the  so-called  artificial  viscosity,  explicit  or  implicit, 
does  not  exist  in  these  schemes. 

We  are  slightly  less  fortunate  with  phase  preservation. 
Exact  solutions  of  the  original  differential  equation  (1) 
should  advect  all  wave  components  through  the  discrete  mesh 
at  precisely  the  speed  u.  Since  we  were  not  able  to  eliminate 
all  the  dispersion  terms  (i.e.,  the  odd  derivatives  of  $ in 
the  Taylor  series)  in  the  BET  procedure,  a certain  amount  of 
numerical  dispersion  or  phase  error  will  occur.  That  is,  for 
a given  size  6x,  Fourier  component  solutions  of  different 
wavelengths  will  be  advected  at  different  speeds.  One  of  the 
most  important  aspects  of  BET  is  the  systematic  reduction  of 
phase  errors  by  eliminating  the  dominant,  lower-order  odd 
derivatives  of  4,  such  as  and  <i>  in  the  Al,  A2,  and  A3 

ya  aAa 

schemes.  In  the  A4  scheme,  the  4/^  term  is  also  eliminated. 
Thus  it  is  expected  to  produce  even  less  phase  error  than  these 
other  schemes. 
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The  analysis  of  phase  error  based  on  Fourier  wave  com- 
ponents is  straightforward,  but  rather  lengthy.  Instead  of 
giving  the  detailed  analyses,  useful  information  can  be 
obtained  by  employing  the  various  schemes  in  actual  computa- 
tions. A standard  test  problem  is  chosen  as  follows:  At  t = 0 

(NCYC  = 0,  where  NCYC  is  the  time  level  cycle  number),  the 
initial  distribution  of  <{>  is  the  Gaussian  function 

<P  = exp{ (In  2)(x  - xQ)2/r2}  (18) 

where  r = 5 and  xq  = 30.  The  computational  domain  is 
0 - x i 60 . To  save  computer  time  in  long  calculations, 
periodic  boundary  conditions  are  imposed. 

To  determine  the  quality  of  a numerical  solution  at  a 
given  point  in  time,  we  simply  compare  it  with  the  exact 
solution  for  the  same  instant.  The  exact  solution  is  identical 
with  the  initial  profile  (Eq.  (18)),  except  with  the  center  xQ 
being  transported  at  the  constant  velocity  u.  In  the  series  of 
tests  conducted  (Figs.  4-7),  the  following  parameters  are 
used 

u ■ 1.0 
6x  = 1.0 

° “ 1ST  " 0,75  * 

The  results  of  Al,  A2,  and  A3  are  compared  with  the 
exact  solution  at  NCYC  - 10,000  in  Fig.  4.  It  is  seen  that 
both  the  Al  and  A2  solutions  are  very  close  to  the  exact 
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solution,  while  the  A3  solution  is  slightly  distorted.  In 
Fig.  5,  the  three  schemes  are  compared  at  NCYC  = 20,000.  Now 
the  deviation  from  the  exact  solution  is  more  visible  for  all 
three  schemes,  with  the  A1  and  A2  results  showing  a small 
lagging  phase  error,  and  the  A3  result  definitely  shows  a 
leading  phase  error.  These  errors  become  quite  pronounced 
at  NCYC  = 50,000  (Fig.  6).  At  this  time,  notice  that  the 
maximum  heights  of  the  numerical  solutions  are  about  4%  to 
6%  lower  than  the  exact  solution.  This  height  reduction  is 
not  a consequence  of  any  numerical  damping,  but  is  rather 
due  to  the  dispersion  of  the  short  wavelength  components  away 
from  the  center  of  the  Gaussian  distribution. 

As  anticipated  earlier,  very  accurate  solutions  can  be 
obtained  by  using  the  A4  scheme.  In  Fig.  7,  the  A4  solutions 
at  NCYC  = 50,000  and  100,000  are  compared,  and  they  are  in- 
distinguishable from  the  exact  solution. 

It  should  be  mentioned  that,  while  A1  and  A4  are  effec- 
tively explicit,  they  become  implicit  when  used  with 
periodic  boundary  conditions.  From  experience,  an  economical 
solution  procedure  is  to  use  A2  as  a predictor  to  find  pro- 
visional  values  for  <frj  , and  then  use  A4  twice  as  correctors. 
The  results  are  very  satisfactory,  as  shown  in  Fig.  7. 

The  Al,  A2,  and  A4  schemes  possess  a very  interesting 
property:  For  a * 0,  1/2,  and  1,  all  the  truncation  errors 
vanich  identically,  as  one  can  readily  verify  by  examining 
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Eqs . (10)  and  (15).  For  those  values  of  a,  these  schemes  are 
exact  and  they  preserve  absolutely  whatever  the  initial  dis- 
tribution of  <f>  is  in  the  domain,  transporting  it  precisely 
at  the  velocity  u.  The  case  a = 0 is  a trivial  one,  while 
for  u = 1 the  mesh  values  of  <p  are  simply  passed  from  point 
to  point.  These  limiting  properties  for  a = 0 and  a = 1 are 
shared  by  many  existing  finite  difference  schemes.  However, 
that  Al,  A2,  and  A4  are  exact  for  a = 1/2  came  as  a little 
surprise.  This  property  may  be  quite  useful,  since  in  most 
applications  6t  may  be  chosen  such  that  a is  in  the  neighbor- 
hood of  1/2. 

Another  interesting  finding  is  that  in  both  Al  and  A2 , 
the  phase  error  is  of  the  leading  type  for  0 £ a < £,  and  it 
becomes  lagging  for  £ £ a £ 1.  Figure  8 illustrates  this 
behavior  with  the  A2  results  at  NCYC  - 50,000,  for  a = 0.25, 
0.5,  and  0.75.  The  curve  for  a = 0.5  is  identical  with  the 
exact  solution. 

Finally,  the  various  schemes  were  put  to  a severe  test 
by  using  a sharp-peaked  triangular  profile,  instead  of  the 
smooth  Gaussian  distribution,  as  the  initial  condition  for  $. 
In  Fig.  9,  the  A2  results  are  compared  with  the  exact  solution 
at  NCYC  = 1400.  For  a = 0.75,  the  lagging  phase  error  typical 
of  A2  is  apparent,  while  for  a » 0.5  the  numerical  solution  is 
identical  with  the  exact  one  (even  the  triangular  peak  is 
preserved) . 
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IV.  BEYOND  THE  COURANT  NUMBER  LIMIT 


To  design  an  accurate  scheme  for  computing  flows  at 
high  Reynolds  number  (see  Section  V for  extensions  to  include 
diffusion),  the  first  requirement,  of  course,  is  that  the 
advection  scheme  itself  does  not  produce  errors  that  over- 
whelm the  effect  of  the  real,  physical  diffusion.  In  this 
regard,  the  Al,  A2 , and  A3  schemes  are  good  candidates  for 
such  applications.  At  very  large  Reynolds  numbers,  more 
accurate  schemes,  such  as  A4 , may  be  required. 

However,  even  a perfect  advection  scheme  is  not  suffi- 
cient for  actual  computation  of  high  Reynolds  number  flows,  if 
the  applicability  of  the  scheme  is  subject  to  the  Courant 
condition,  Eq . (11),  which  limits  the  maximum  size  of  fit  to  be 

«t  * ^ . (19) 

In  explicit  calculations,  the  stability  restriction  on  fit  for 
diffusion  {1}  is 


6t  ' I (20) 

where  v is  the  diffusion  coefficient.  At  large  Reynolds  num- 
bers, Eq . (19)  is  far  more  restrictive  than  Eq.  (20).  For 

— h 2 

example,  consider  a flow  with  u = 1 ft/sec  and  v = 10  ft  /sec. 
Assuming  the  characteristic  length  of  the  flow  to  be  D = 10  ft, 
and  fix  - 1 ft,  the  Reynolds  number  of  the  flow  is  then 


Re 


f.  = 106 
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which  is  a very  common  flow  regime.  By  Eq . (20),  one  would 
be  allowed  to  use  fit  up  to  5 * 104  sec;  however,  Eq . (19) 
limits  fit  to  be  less  than  1 sec.  The  consequence  is  that, 
even  using  an  extremely  accurate  advection  scheme,  the  Courant 
condition  would  require  millions  of  time  steps  to  obtain  any 
meaningful  results  for  the  example  above. 

All  the  schemes  presented  in  this  paper  can  be  gener- 
alized to  arbitrary  Courant  number.  We  shall  illustrate  the 
principle  by  considering  the  A1  scheme.  Referring  to  Fig.  10, 
let  s be  a positive  integer  (including  zero)  such  that 

s a a < s + 1 . (21) 

Also,  let 


so  that 


a'  = a - 8 


0 * a'  < 1 

For  any  positive  values  of  a,  by  Eq.  (2)  4/:  is  equal  to 

which  is  located  at  a distance  afix  to  the  left  of  4> j . Com- 
paring the  six  quantities  <t>^,  <J>B,  $c,  4>D,  4»E,  and  4>F  as  shown 
in  Fig.  10  with  those  in  Fig.  1,  one  can  easily  write  down, 
in  analogy  with  Eq.  (10),  the  generalized  A1  scheme  as 


Generalized  A1  Scheme 


.n+i  . n-1 

+ a 


'Ln+1 

ih-i 


J 


(22) 


where 
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,>  „ (1  - a'Hl  - 2a'] 


‘1  (1  + a')(  1 + 2a') 


- _ 2(1  - 2a'] 


Similarly,  the  generalized  forms  for  the  other  schemes  are 


Generalized  A2  Schemes 

An+1  .n-1  . . *f.n  .n  ) . „ f .n  .n  \ 

♦j  *j-2s-l  + bll*J-s-2  * *j-s+lj  + b2(*j-s  " ♦j-s-lJ 


i «'(1 


- a')(l  - 2o') 


bj  - j (i  - 2o')(l  + 0(2  - o') 


Generalized  A3  Scheme 


An+1  .n  . An+11  + „'Ln 

*J  ^ J —8—1  + cl[*i-s  “ ♦j-lj  + c2 [^j-s-2  " 


2 - o' 


c2  “ 


Generalized  A4  Scheme 


♦"+1  ‘ + 4"!i  - ♦3-J.)  ♦ i - 2) 

♦ 4;-.  - 


(25) 
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where 


d 


1 


I 


2 - a 
2 + g 


1 - 2a') 
1 + 2a; 


(Q2  (1  - 2a-) 

3(2  + a') 

(2  - g')(l  - 2a') 


As  in  their  original  forms,  the  generalized  Al,  A2, 
and  A4  schemes  are  exact  for  a'  = 0,  i,  and  1.  Calculations 
have  been  performed  with  g = 100.75  and  1000.75  for  the 
problem  discussed  in  connection  with  Eq . (18).  Each  computa- 
tion was  continued  for  at  least  10,000  time  steps,  and  the 
quality  of  the  results  is  consistent  with  those  mentioned  in 
Figs.  4-7.  Therefore,  unlike  many  unconditionally  stable 
schemes,  the  use  of  large  6t  in  Eqs . (22)-(25)  does  not  impair 
the  accuracy  of  the  solutions.  In  fact,  the  only  parameter 
that  affects  the  accuracy  (i.e.,  phase  preservation)  is  a', 
the  size  of  g being  irrelevant. 
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V.  ADVECTI ON -DIFFUSION  SCHEMES 


The  model  equation  for  one-dimensional  advection- 
dif fusion  process  is 

S ♦ « & - * § <26> 

where  v is  a constant  kinematic  viscosity  or  diffusion  coeffi- 
cient. The  numerical  solution  procedure  will  be  illustrated 
by  extending  the  basic  A1  scheme  to  include  diffusion. 

The  resulting  discrete  representation  of  Eq . (26)  is  called 
the  ADI  scheme . 

5.1  The  ADI  Scheme 

Consider  Fig.  11,  which  may  be  compared  with  Fig.  1 

for  analogy.  Unlike  in  Fig.  1,  is  not  equal  to  4>j+1 

which  is  the  value  of  + located  at  a distance  afix  to  the  left 

of  Although  the  fluid  particle  originally  situated  at 

the  location  of  does  move  to  the  position  of  4>”+1  at  the 

end  of  the  time  increment  fit,  the  diffusion  process  has 

changed  the  $ value  associated  with  this  particle.  From 

this  Lagrangian  interpretation,  one  can  relate  and 

by  an  explicit,  forward-time,  central-space,  finite  difference 

2 2 

scheme  for  the  diffusion  equation  D$/Dt  = v( 3 $/3x  ),  where 
D/Dt  is  the  Lagrangian  derivative.  That  is, 

♦r1  ■ *s+i  ♦ »(♦$  - 2 ♦r1  + ♦3-i)  (2?) 
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where 


For  stability,  it  is  required  that  0 £ 3 £ 1/2  { 1 } . 

Similarly,  because  of  diffusion,  ♦j-1  is  not  equal  to 
Rather,  they  are  related  by  the  explicit  formula 

♦yl  - *vx  * »(♦;;!  - * ♦r1  * *v-i]  <28> 

In  close  analogy  with  Eq . (10),  the  six  quantities 

♦j-1’  ♦j-l»  ♦jll*  *y  »nd  ♦5"1  (Tig-  H>  sre  expanded 

in  Taylor  series  about  to  give 


+ a. 


n+1 

j-1 


+ 


a2 


(29) 


where  aj  and  a2  are  the  same  as  given  in  Eq . (10). 

The  solution  procedure  in  the  ADI  scheme  consists  of 
three  steps: 

~n-l 

1.  Calculate  , using  Eq.  (28),  for  the  entire 
space  domain. 

2.  By  Eq.  (29),  calculate  4»“+1  throughout  the  domain. 

3.  Use  Eq.  (27)  to  obtain  *"+1. 

The  procedure  above  is  explicit  if  in  Step  2 is  computed 

in  the  direction  of  increasing  j . 

The  ADI  scheme  can  also  be  generalized  to  allow  arbi- 
trary Courant  number,  and  the  resulting  procedure  is 
1.  Calculate  throughout  the  domain  by 


rn+l 

♦j 


in_l 

♦j-2s-l 


♦S3.)  + »a(*J-s  - ♦j-s-l)  <30> 
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The  "half-strength”  radius  r is  the  distance  from  x = xQ  such 
that  <J>(xo+r,0)  = 1/2.  For  this  test  we  used  r = 5 ft,  thus 
the  characteristic  length  D = 2r  = 10  ft.  The  Reynolds  number 
is  then 

Ite  = — = 10^  . 

V 


is 


The  exact  solution  which  satisfies  Eqs . (26)  and  (34) 


4>(x,t) 


(x  - xQ  - ut)2 
4v(t  + tQ) 


(36) 


A mesh  of  100  6x  (fix  = 1 ft)  in  length  was  used,  the 
point  x = xQ  being  at  the  center  of  this  domain  at  t =0. 
Periodic  boundary  conditions  were  imposed  to  economize  the 
computation.  For  diffusion  stability,  Eq . (20)  requires  that 
fit  £ 5 x lo3  sec.  We  used  fit  = 1000.75  sec,  which  makes  the 
Courant  number  a = 1000.75.  The  AD3  scheme,  an  extension  of 
A3,  was  used  to  start  the  first  time  step.  Afterwards,  the 
ADI  scheme  was  employed. 

In  Fig.  12,  the  numerical  solutions  are  compared  with 
the  exact  solution,  Eq.  (36),  for  NCYC  = 100  (t  = 100075  sec) 
and  NCYC  * 360  (t  = 360270  sec).  The  agreement  is  excellent. 
As  mentioned  earlier,  if  one  were  restricted  by  the  usual 
Courant  condition  0 < a 1 1,  the  maximum  allowable  fit  would 
be  less  than  1 sec.  It  would  then  take  more  than  100075  time 
steps  to  obtain  the  solution  at  t = 100075  sec,  which  would 
make  the  cost  of  computation  a rather  serious  matter. 
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5 . 3 Boundary  Conditions 

In  the  sample  calculations  mentioned  so  far,  periodic 
boundary  conditions  have  been  employed.  There  are,  however, 
numerous  advection-dif fusion  flow  situations  where  other  types 
of  boundary  conditions  arise.  A frequently  encountered 
boundary  condition  is  the  specification  of  <p  at  the  upstream 
and  downstream  ends  of  the  domain.  We  shall  limit  our  dis- 
cussion here  to  boundary  conditions  needed  in  the  basic  ADI, 
AD2,  and  AD3  schemes  (i.e.,  the  Al,  A2 , and  A3  schemes  gen- 
eralized to  include  diffusion,  but  only  for  0 £ a £ 1) . 

First,  consider  the  downstream  boundary  (at  j = R in 
Fig.  13).  The  boundary  quantities  needed  for  computation  in 
the  interior  are  and  <j>R.  Since  4>R  is  specified,  we  only 

have  to  find  an  expression  for  $R+1.  This  is  done  by  making 
balanced  expansions  of  <j>R  $R  » and  <J>R  about  <t>Q 

(Fig.  13).  The  result  is 


c1  ■ *ti + (Hi)  K - c;) 


(37) 


The  situation  at  the  upstream  boundary  (j  = L in  Fig.  13) 
is  more  complicated.  The  boundary  quantities  needed  for  com- 
puting at  the  interior  point  j = L+l  are 

and  which  is  known.  The  first  step  is  to  find  <f>£_ ^ . We 
write  Eq . (26)  in  a forward-time,  central-space  difference 


about 
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.n+1 


An  Ln  An 

r-*k  ♦ uN1 ' 


6t 


26x 


'L+l  " 20L  + *L-H 
6x2 


which  can  be  rearranged  as 


PL-1 


(a  - 2B)  <J>£+1  + 4B  4>l 


a + 2B 


.n+1 


- *i 


(38) 


The  second  step  is  to  calculate  1 by  Eq • (28)-  To  do  this, 
we  need  4>l_i  which  can  be  obtained  from  Eq . (38)  by  reducing 
all  superscripts  by  one.  The  third  step  is  to  find  $£+1. 

Using  $L+1'  ^L’  and  $l_1  (Fig.  13),  we  again  form  a 

forward-time , central-space,  difference  approximation  to 
Eq.  (26)  as 


An+1 

,n  /j.n-1  i-n+l\ 

*L  + /* L - *L  ) 

- vl 

fe-1  - k * clNi 

6t 

l 2a  fix  J 

v! 

( a2  6x2  / 

By  rearranging,  we  obtain 


,„+l  <“2  - 28)  if1  * 46  + 2a‘ 

>L  ^~n~e 


.n+1 


- <t> 


n 


(39) 


The  preceding  steps  provide  all  the  boundary  conditions  needed 
for  calculating  + 1 at  j = L+l. 

5 .4  Steady-State  Solution  and  the  Cell  Reynolds  Number  Limit 
To  validate  further  the  advection-dif fusion  schemes 
presented  here,  we  have  applied  them  to  the  steady -state 
problem 


27 


u|i=  vili 


9x 


9x 


(40) 


subject  to  the  boundary  conditions 

<j>  = 0 at  x = 0 

4>  = 1 at  x = 1 

The  exact  solution  to  this  problem  is  {2} 

1 _ e(ux)/v 

♦<*>  - T-V"  ' 


(41) 


The  ADI  scheme  with  the  boundary  conditions  treated  in 
the  preceding  section  was  employed.  A set  of  finite  difference 
equations  for  the  steady  state  was  obtained  from  ADI  by  setting 
= ^j”^  = for  values  of  j.  Following  Roache  (2), 
we  used  fix  = 0.1  and  solved  the  steady -state  equations 
by  direct  Gauss  elimination.  The  numerical 
results  are  compared  with  the  exact  solution  in  Fig.  14,  for 
v/u  = ®,  1.0,  0.25,  0.10,  and  0.05.  For  v/u  less  than  0.10, 
the  numerical  solution  becomes  less  accurate  as  v/u  decreases, 
because  of  inadequate  resolution  by  the  mesh.  The  accuracy 
can  be  improved  by  increasing  the  resolution.  In  Fig.  15, 
numerical  solutions  using  fix  ■ 0.01  are  compared  with  the  exact 
solutions.  The  agreement  is  excellent. 

A difficulty  suffered  by  many  existing  advection- 
dif fusion  schemes  is  the  "cell  Reynolds  number  limit,"  which 
has  been  described  in  detail  by  Roache  {1,  2,  3).  The  cell 


28 


r 


Reynolds  number  (Rc)  is  the  local  Reynolds  number  based  on  6x 
as  a characteristic  length,  viz. 


R = 

C V 


(42) 


In  many  schemes,  e.g.,  central  difference  for  3<(>/3x  and 
2 2 

3 <f>/3x  in  Eq . (40),  the  steady-state  numerical  solutions  are 
qualitatively  similar  to  the  exact  solutions  for  R 1 2.  But 
when  R > 2,  the  discrete  solutions  are  oscillatory  and  diverge 
badly  from  the  exact  solutions. 

Using  any  one  of  the  advection-dif fusion  schemes 
presented  in  this  paper  (e.g.,  ADI),  solutions  have  been  obtained 

O 

for  R > 2 and  up  to  R = 10  . Table  I shows  the  numerical 

C G 

results  for  those  cases  using  fix  = 0.1.  The  discrete  <(>  values 

5 

shown  in  the  table  have  been  multiplied  by  a factor  of  10  for 
easy  reading.  It  is  seen  that  as  R increases  (or  as  v/u 

G 

decreases  for  a fixed  fix),  $ becomes  vanishingly  small  over  a 
greater  portion  of  the  domain  and  it  rises  more  sharply  near 
x 3 1 to  meet  the  downstream  boundary  condition  4>(  1 ) = 1.  At 

Q 

R « 10  , a literally  discontinuous  solution  is  obtained,  which 
c 

is  in  agreement  with  the  exact  solution. 

Since  fix  and  fit  appear  in  the  steady-state  difference 
equations  derived  from  ADI,  a question  may  be  raised  with  regard 
to  the  uniqueness  of  the  numerical  steady-state  solution.  That 
is,  for  a given  pair  of  u and  v (u  f 0 and  v f 0),  only  one 
exact  solution  exists,  according  to  Eq.  (41).  But  in  a 


numerical  model,  the  solution  may  depend  on  the  choice  of  the 
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nonphysical  parameters  6x  and  fit.  Once  a suitable  fix  has  been 

selected,  then  for  the  given  u and  v the  only  free  parameter 

is  fit.  Thus  the  question  boils  down  to  whether  the  steady-state 

numerical  solution  is  sensitive  to  variation  in  fit.  For  the 

various  calculations  shown  in  Figs.  14,  15,  and  Table  I, 

identical  results  (i.e.,  four  significant  figures)  were  obtained 
—5 

for  10  < a ^ 0.99.  For  practical  purposes,  therefore,  the 

present  steady-state  solutions  are  independent  of  fit,  as  they 


should  be . 


VI.  CONCLUSIONS 


The  various  numerical  experiments  conducted  so  far 
indicate  that  the  advection  and  advection-di f fusion  schemes 
presented  in  this  paper  are  highly  accurate.  In  particular, 
the  advection  schemes  are  all  neutrally  stable,  free  of  any 
artificial  viscosity  (positive  or  negative).  The  A1  and  A2 
schemes  produce  little  phase  errors,  and  accurate  results  have 
been  obtained  up  to  20,000  time  steps  (at  a = 0.75)  by  using 
either  method.  The  A3  scheme  is  fully  implicit  and  more 
costly  to  use,  but  it  is  needed  only  for  the  first  time 
step  as  a starter  for  the  other  schemes.  The  most  accurate 
method  considered  so  far  is  the  A4  scheme,  which,  after 
100,000  time  steps  (at  a = 0.75),  yields  results  that  are  hardly 
distinguishable  from  exact  solutions.  More  importantly,  the 
success  of  these  schemes  demonstrates  that  for  some  differen- 
tial equations,  such  as  Eq.  (1),  evaluation  of  the  derivatives 
is  not  the  only  way  to  obtain  finite  difference  representations 
of  the  original  differential  equation. 

With  the  removal  of  Courant  number  limit,  the  present 
schemes  can  be  used  to  compute  flows  at  very  high  Reynolds 
number,  at  least  for  the  model  equation  considered.  Further- 
more, with  sufficient  spatial  resolution  the  present  advection- 
diffusion  schemes  produced  steady-state  solutions  that  not  only 
agree  with  exact  solutions,  but  are  completely  free  from  the 
cell  Reynolds  number  limit. 
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While  it  is  true  that  Eqs.  (1)  and  (26)  represent  fairly 


simple  systems  (i.e.,  linear,  with  constant  coefficients,  and 
only  for  one  space  dimension),  they  are  nevertheless  important 
for  testing  the  validity  of  any  proposed  advection  or  advection- 
dif fusion  schemes.  By  the  general  principle  described  in  this 
paper,  the  author  believes  that  accurate  schemes  can  be  derived 
for  multi-dimensional  problems,  flows  with  nonuniform  velocity 
and  viscosity,  and  mesh  with  variable  spacings.  These  ex- 
tensions are  currently  under  study,  and  the  findings  will  be 
reported  at  a future  date. 


I 
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Fig.  1.  Definition  of  the  Variables  Used  in  A1  Scheme. 


Comparison  of  Schemes  at  NCYC  **  10,000,  using  a = 0.75 


of  Schemes  at  NCYC  ■ 20,000,  using  a = 0.75 


Fig.  12.  Comparison  of  ADI  Results  (dots)  with  Exact  Solution 


Fig.  13.  Definition  Sketch  for  Boundary  Conditions. 


Exact  Solution  (solid  lines)  . The  spatial  resolution 
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